const char* dgemm_desc = "Naive, three-loop dgemm.";

/*#include <memory.h> */
#include <stdlib.h>
#include <stdio.h>
//#include <xmmintrin.h>
//#include <emmintrin.h>
#include <x86intrin.h> // includes all SSE intrinsics!
#include <string.h>
#include <malloc.h>
/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. */    
void square_dgemm (int n, double* restrict A, double* restrict B, double* restrict C)
{


    /* Reorganize A to be row-major */
    /* TODO: How to do this efficiently? */
    double* Atemp = (double *)malloc(sizeof(double)*n*n);
    //void ** restrict Atemp;
    //posix_memalign(Atemp,128,sizeof(double)*n*n);
    (double *) Atemp;
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            Atemp[i*n+j] = A[i+n*j];
        }
    }
    /* Also copying B so that it's aligned */
    double* restrict Btemp = (double *)malloc(sizeof(double)*n*n);
    memcpy(Btemp,B,sizeof(double)*n*n);

  /* For each column j of B */
  for (int j = 0; j < n; ++j) 
    /* For each row i of Atemp */
    for (int i = 0; i < n; ++i)
    {
      /* Compute C(i,j) */
      double cij = C[i+j*n];
      double* restrict Arowi = Atemp+(i*n);
      double* restrict Bcolj = Btemp+(j*n);
    int k=0;
    for ( ; k<n-8; k+=8) {
	    // Try loop unrolling
            cij += Arowi[k] * Bcolj[k];
	    cij += Arowi[k+1]*Bcolj[k+1];
      	    cij += Arowi[k+2]*Bcolj[k+2];
	    cij += Arowi[k+3]*Bcolj[k+3];
	    cij += Arowi[k+4]*Bcolj[k+4];
	    cij += Arowi[k+5]*Bcolj[k+5];
	    cij += Arowi[k+6]*Bcolj[k+6];
	    cij += Arowi[k+7]*Bcolj[k+7];
	}
	if (k!=n){
        for (;k<n;k++){
            cij += Arowi[k]*Bcolj[k];
        }
    }

      C[i+j*n] = cij;
    }

    //It's always a good idea to free allocated memory
    free(Atemp);
    free(Btemp);
}
